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Abstract. With the aim to get an insight in the origin of differences in the ear- 
lier reported calculation results for KNbC>3 and to test the recently proposed "NFP" 
implementation of the full-potential linear muffin-tin orbital (FP-LMTO) method by 
M. Methfessel and M. van Schilfgaarde, we perform a comparative study of the ferro- 
electric instability in KNbC>3 by FP-LMTO and full-potential linear augmented plane- 
wave (LAPW) method. It is shown that a high precision in the description of the 
charge density variations over the interstitial region in perovskite materials is essen- 
tial; the technical limitations of the accuracy of charge- density description apparently 
accounted for previously reported slight disagreement with the LAPW results. With 
more accurate description of the charge density by sufficiently fine real-space grid, the 
results obtained by both methods became almost identical. 

In order to extract additional information (beyond the harmonic approximation) 
from the total energy fit obtainable in total-energy calculations, a scheme is proposed 
to solve the multidimensional vibrational Schrodinger equation in the model of non- 
interacting anharmonic oscillators via the expansion in hyperspherical harmonics. Pre- 
liminary results are given for the t\ u vibrational modes in cubic KNb03. 



INTRODUCTION 

KNbOa is one of benchmark systems for ab initio analysis of ferroelectric per- 
ovskites. It has been extensively studied by the whole spectrum of numerical 
methods from an apparently ultimately accurate full-potential linear augmented 
plane waves (LAPW) [1-5] through pseudopotentials [6] and tight-binding ab ini- 
tio schemes [7] to semiempirical and model-based techniques [8] . 

Modern state-of-art simulations of ferroelectric correlations, lattice dynamics and 
phase transitions are dependent on reliable and accurate description of the total 
energy as function of displacements and strain variables. The full-potential linear 
muffin-t in-orbit als (FP-LMTO) method [9] was proven able to provide a reasonable 
balance of accuracy and low computational effort even when applied to supercells of 
up to 40 atoms [10-13]. However, in sensitive benchmark calculations of T phonons, 



due to certain technical limitations, FP-LMTO provided apparently lower accuracy 
for phonon eigenvectors, as compared to FP-LAPW [11,3,5]. A new version of the 
FP-LMTO code by Methfessel and van Schilfgaarde [14], in contrast to the pre- 
vious one used in our earlier calculations [9], is unsensitive to sphere packing and 
uses much more efficient, albeit more mathematically involved, basis of augmented 
"smooth Hankel functions" that enables one to drastically reduce the size of di- 
agonalization problem without loss of accuracy. We compare the results obtained 
with the new FP-LMTO and with the WIEN97 implementation of the FP-LAPW 
method [15], concentrating on the accuracy and performance. 

In order to get an additional insight into lattice dynamics of KNbOs, we present 
vibrational frequencies as calculated quantum-mechanically in the assumption of 
uncoupled multidimensional anharmonic oscillators, based on the total energy data 
obtained from first-principles calculations, and give a preliminary estimate of the 
lowest vibration frequency within this approach. 

COMPARISON OF CALCULATION METHODS 

The linear augmented plane-wave (LAPW) and the linear muffin-tin orbital 
(LMTO) methods are closely related and originate in the same work by Ander- 
sen [16]. In the modern stage at their development, they share the advantage of 
being all-electron methods (in contrast to pseudopotential ones and to those de- 
pending on the frozen-core approximation). As compared to tight-binding schemes 
with fixed (e.g. Gaussian-type) bases, the basis in LAPW and LMTO is optimized 
in the course of iterations, as the heads of wavefunctions are recalculated inside the 
(arbitrarily) chosen muffin-tin spheres from one iteration to another. The tails of 
basis functions that span the interstitial region between the spheres are constructed 
by an augmentation procedure which matches the numerical solutions inside the 
MT-spheres either to the plane waves of different k (in LAPW) or to the spherical 
Hankel functions of different energy (in LMTO). This difference accounts for rel- 
ative advantages and disadvantages of these two calculation schemes. In LAPW, 
the number of augmented plane waves with different k can be safely saturated, 
until the desirable convergence of results with respect to the completeness of basis 
is achieved. This however has its price in terms of computational effort. 

The advantage of LMTO lies in the fact that its basis functions may be tuned to 
resemble true atom-centered wavefunctions in a crystal. This means, in the ideal 
case, quite efficient and compact basis set and hence computational speed and the 
ability to treat larger systems than is possible with LAPW, given certain amount 
of computer resources. The weak point of LMTO is that the saturation of the basis 
set is not as straightforward as in the LAPW, because adding more atom-centered 
tail functions of certain angular symmetry needs care to prevent linear dependences 
within the basis set. In order to overcome this problem and at the same time to 
maintain the LMTO as, ideally, the 'minimal basis' method of competitive accuracy, 
one can try to experiment with more sophisticated but hopefully more efficient 



envelope functions. As some examples of proposed alternatives to conventional 
spherical Hankel functions one can single out the tight-binding LMTO method [17] 
or the "exact muffin-tin orbital theory" development [18], both of them being not 
yet implemented, to our knowledge, in workable full-potential total-energy codes. 

Recently, yet another extension within the general LMTO formalism has been 
proposed and implemented by M. Methfessel and M. van Schilfgaarde and referred 
to as the "NFP" LMTO code [14]. While in many aspects a development of the 
earlier described [9] and widely used FP-LMTO formalism, the NFP algorithm 
incorporates an essential new element, that is, using the "smooth" Hankel functions 
rather than the "standard" ones for the augmentation of numerical radial solutions 
in the interstitial region. Whereas the standard spherical Hankel function satisfies 
the differential equation 

(A + e)ho(r) = -47r5(r) 

for / = 0, the equation for the smooth function contains the smeared 5 function 
go(r) = C exp(— r 2 /i?g m ) as a source term 

(A + e)~h (r) = -A7ig {r) , 

and the smooth functions of higher orders are generated by applying a differential 
operator y(— V), defined by y(r) = r l Y L , to the function of the 0th order. These 
new envelope functions can be tuned by a proper choice of the smoothing parameter 
i?sm to imitate the actual shape of the wavefunction in crystal. 

Cubic perovskite ferroelectrics provide an excellent benchmark system for the 
high-precision total-energy calculation scheme. Whereas one typically cannot pin- 
point any essential disagreement in the band dispersions calculated by different 
full-potential schemes, the energy differences on the ~1 mRy scale related to fer- 
roelectric instability, soft mode phonon frequencies and eigenvectors in these ma- 
terials may be quite differently estimated by different computation schemes, result 
in qualitatively different predictions and thus dramatize the competition between 
different numerical approaches. For KNbOs, a number of calculations has been 
already done using different methods. Full-potential LAPW calculations were per- 
formed by Singh [1-3] and Krakauer et al [4,5], LMTO calculations by us [10-13]. 
In spite of the overall agreement (the ability of both methods to account for the fer- 
roelectric instability, consistent results for the T-frozen phonon frequencies), some 
disagreements prevailed in the description of ferroelectric instability, as well as in 
the estimation of the soft mode phonon eigenvector [11,3,5]. 

What makes perovskites generally (and KNb0 3 as an example) a hard test for 
any computational scheme relying on a site-centered basis, like LMTO, is their 
relatively loose structural packing (if one thinks in terms of nominal ionic radii). 
The electron density is unevenly distributed between compact Nb06 octahedra 
and intermediate large cavities which host small K ions. In our earlier LMTO 
calculations using the "old" LMTO code by Methfessel [9], good sphere packing 
was essential for accurate integration over the interstitial region, but could not be 
guaranteed in a completely satisfactory way (see, e.g. the discussion in Ref. [13]). 



TABLE 1. NFP-LMTO calculation setup (k 2 in Ry, i? MT ,SM in a.u.) 
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In the present implementation of LMTO, the aspect of good sphere packing is 
no more sensitive, therefore we tried to figure out finally whether the mentioned 
disagreement with the LAPW results was due to problems of technical inefficiency 
(unadequate integration scheme etc.) or has to do with some basic limitations of 
the LMTO formalism. 

In our present benchmark calculations, we used the implementation of the full- 
potential LAPW method knows as the WIEN97 code [15]. Sphere sizes were chosen 
as shown in the Table 1 (the same for LAPW and LMTO). The k-space integration 
was performed in an identical way in both schemes, using the sampling on a mesh of 
18 inequivalent points, corresponding to 6x6x6 divisions of the full Brillouin zone. 
This mesh was found to be sufficiently dense for the estimations of ferroelectric 
instability, based on previous experience [2] . The exchange-correlation was treated 
either in the local density approximation (LDA) or in the generalized gradient 
approximation (GGA). 

In the construction of the LMTO setup, the usual procedure is to optimize the 
possibly minimal basis, using the freedom in the choice of basis parameters, and 
then to expand the basis in order to ensure its sufficient completeness. In the 
NFT formalism, the quality of the basis depends on both energies of the smooth 
Hankel functions and their smoothing radii, the latter apparently having more 
pronounced effect. In contrast to earlier FP-LMTO scheme [9] which favored the 
basis functions with at least three tail energies per orbital for sufficient accuracy, the 
NFP provides a reasonable description of the valence band states with augmenting 
just one smooth Hankel function to each; the setup is then refined by adding some 
other tail functions. The calculation setup we used in the present LMTO calculation 
is shown in Table 1. It includes 70 basis functions, that can be either somehow more 
extended, or reduced, to get the desirable compromise between the accuracy and 
performance. For comparison, the basis size for a LAPW calculation of comparable 
accuracy should include at least ~800 augmented plane waves. 

The present version of NFP incorporates only the LDA treatment of the ex- 
change/correlation. Another technical drawback is the impossibility to treat the 
states with different principal quantum number and the same orbital quantum 
number within the valence band. For KNb0 3 , we had to include Nb4p states and 
neglect Nb5p in the valence panel. This seems to be an acceptable compromise, 
however. More discussion related to the LMTO setup may be found in Ref. [10]. 



CALCULATION RESULTS AND DISCUSSION 



Fig. 1 shows the energy/volume curves as calculated by LAPW in the LDA and 
in the GGA, and by LMTO in the LDA. First two curves essentially reproduce 
previous results by Singh [2] aimed at the comparison of LDA and GGA. The 
energy/volume curve generated now with LMTO (crosses in Fig. 1) practically 
coincides with that obtained with the LAPW. Absolute energy values lie by ~0.9 
Ry lower with LAPW, understandably due to more complete basis. 

For the study of ferroelectric instability, we concentrated on the displacement 
pattern compatible with the t± u TO phonon modes, i.e. the ^-displacements of 
K(0 0), Nb(iii) and 0„ (||0) with respect to two equivalent : (0±±) and (0±±) 
atoms. In Fig. 2, the energy differences are shown as function of the displacement 
pattern which roughly corresponds to that in the soft mode, ultimately resulting 
in the equilibrium structure of tetragonal ferroelectric phase: Nb is displaced twice 
farther as K relatively to the oxygen cage. We found the calculated energy differ- 
ences to be extremely sensitive to the quality of the charge density representation 
in the unit cell. In both methods we used, this expansion is done by the fast Fourier 
transformation; in WIEN97, the magnitude of the largest reciprocal-space vector 
G is specified whereas in the NFP-LMTO the number of divisions N along each 
unit cell edge for a real-space uniform grid has to be explicitly provided. For a per- 
ovskite, both cutoffs need to be relatively large in order to achieve a convergency in 
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FIGURE 1. Total energy difference (with respect to the equilibrium value) depending on the 
lattice constant in KNb0 3 as calculated by LMTO-LDA (crosses), LAPW-LDA (filled circles) 
and LAPW-GGA (open circles). The parabolic fit is shown for each case. Experimental lattice 
constant extrapolated to zero temperature is indicated by a vertical line. 
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FIGURE 2. Total energy difference depending on atomic displacements A Z (K) + 2A z (Nb) as 
calculated by LAPW and LMTO for different values of A z and different cutoffs in the charge 
density expansion. Crosses: G=12.0 in LAPW and N=18 in LMTO; dots: G=U.O in LAPW 
and N=24 in LMTO (essentially converged results). Polynomial fit is a guide to eye. 



TABLE 2. Calculated frequencies and eigenvectors of T-TO phonons of the t\ u symmetry from 
LAPW and LMTO 
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this parameter. As is seen in Fig. 2, the value G = 10 in LAPW overestimates the 
ferroelectric instability whereas the LMTO with = 18 finds yet no trace of the 
instability. At G — 12 and iV = 24, the trends are about comparable, becoming 
even closer for G = 14 and N = 32. 

As an additional test for the proper balance of the energetics of different displace- 
ment patterns, we calculated the T-TO phonons in the cubic phase of KNb0 3 . The 
insufficient accuracy of previous LMTO calculations [11,12] manifested itself as a 
noticeable deviation of the vibration eigenvectors from those obtained by LAPW 



[1,5]. Table 2 shows the results of some preliminary calculations with the NFP 
code (obtained with the 2d-order total energy fit over the results for several com- 
bined displacements) in comparison with the LAPW data. One can see that the 
correct displacement pattern within the soft mode is now restored, and the overall 
agreement with the LAPW eigenvectors is quite satisfactory. 



TREATMENT OF ANHARMONIC VIBRATIONS 

With the total-energy fit generally available from first-principles calculations, 
one may tend to extract some additional information than is possible within the 
harmonic approximation. The treatment of anharmonic effects in crystal is rather 
complicated (see, e.g., Ref. [19] for a review). In principle, the modes of different 
symmetry and related to different q- values couple beyond the harmonic approxima- 
tion. Nevertheless, in the study of ferroelectrics there have been several attempts 
to single out any particular mode, which is believed to be principally associated 
with anharmonic effects, and to solve the vibrational Schrodinger equation related 
to it. This has been done e.g. for LiTaC>3 and LiNb03 by Inbar and Cohen [20] 
and by Bakker et al. [21] (for an empirical potential well, in the latter case) as 
well as for a one-dimensional A2 mode in orthorhombic KNbOs by Postnikov and 
Borstel [12]. This approach was referred to as non- interacting anharmonic oscil- 
lators [20], meaning the oscillators related to a particular T TO-mode in crystal. 
It is assumed that the displacement potential is separable into components with 
different q-dependence. Such separation is less valid for several symmetry coordi- 
nates which mix already in the harmonic approximation, therefore the solution of 
a multidimensional oscillator problem is necessary in this case. A straightforward 
treatment by, e.g., a finite-difference method on a multidimensional grid rapidly 
becomes prohibitive with a number of dimensions (see [12]). Therefore, we propose 
a scheme which uses the expansion in hyperspherical harmonics. This approach is 
known in the calculation of vibration spectra of three-atomic molecules [22], how- 
ever, for an ad hoc constructed system of variables. In the following, we describe 
the formalism for an arbitrary number of symmetry coordinates. 

We start with an arbitrary convenient set of symmetry-adapted displacement 
coordinates (see, e.g., Ref. [23]): S t = J2iBuXi (xi are conventional cartesian dis- 
placements), which form a complete basis within a particular irreducible represen- 
tation, but do not need to be orthonormal. The Schrodinger equation then acquires 
a form: 



^ = £^, (1) 



with the kinetic-energy matrix G tt ' = J2i B t im i 1 B t 'i- Mixed derivatives can further 
be excluded by the following orthogonalizing transformation: 



where X t / t is the t-th eigenvector, corresponding to the eigenvalue A*, of the kinetic- 
energy matrix. In n-dim. space of generalized coordinates Q t , we use spherical 
coordinates (r, ■ ■ ■ , $ p , <p) for p = n — 2, 
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The explicit form of the polynomials of is the following: 
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The hyperspherical harmonics are chosen either as complex functions 
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or in the real form, ~ cosm p for m p > and ~ sinm p for m p < 0. C£(z) are 
Gegenbauer polynomials [24,25], generated by the following recursion: 

C p (z) = 1; C?(z) = 2pz; (n+1) C£ +1 (z) = 2(n+p) z C^(z)-(n+2p-l) C^z) . 

^ _ i (t?i,.,t? p , </?), orthogonal on a unit sphere, are eigenfunctions of the multi- 
dimensional Laplace operator: 

AY - m °( m °+P) y. 

With the potential and the wave function expanded in hyperspherical harmonics 
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the multidimensional Schrodinger equation (1) transforms into a system of 
X^=o x '(2Af + p)(A^+p— l)!/(p! AT!) coupled 1-dimensional equations: 
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The expansion of the potential (provided in a polynomial form by a fit to total- 
energy values) is finite whereas for the wave function a cutoff value iV max has to be 
introduced. 

As a practical example of this approach, we considered the solution of a 3- 
dimensional oscillator problem corresponding to the vibration pattern within the 
t\ u mode in cubic KNb0 3 . For the symmetry coordinates as discussed above, we 
included the 4th power of the Nb displacement into the total energy fit. (For real 
applications, one should of course consider some other degrees of freedom beyond 
the harmonic approximation). The system of coupled equations (2) was solved by 
a finite difference method, with 50 points in the equidistant radial mesh from up to 
r = 5.0 where a boundary condition R(r) = was imposed on radial wavefunctions 
(this scheme may be somehow refined in more precise calculations, incorporating a 
nonuniform mesh). For the maximal degree of polynome iV = 8 in the wave func- 
tion expansion, the energy difference between two lowest oscillator levels practically 
converged to 70 cm -1 . The convergence is more slow for higher levels. 

CONCLUSIONS 

We compare in the present paper the results obtained for the ferroelectric insta- 
bility in KNbOs with two methods, LAPW and LMTO, both of which have been 
applied to this system before but apparently never underwent a thorough com- 
parison with the, as far as possible, identical calculation setup. The result of this 
comparison is that not only the energy/volume curves are identical in the LDA, but 
the description of the ferroelectric instability, involving equilibrium displacements 
of ~0.1 a.u. and energy differences of ~0.5 mRy, is practically identical by both 
schemes, provided the sufficient accuracy in the description of the charge density 
variations over the unit cell is guaranteed. The LAPW method provides under- 
standably lower absolute values of the total energy, but the new formulation of 
LMTO has the advantage of much more compact basis set (about 10 times smaller 
than that of LAPW) and some resources to expand the basis somehow for even 
better controllable accuracy without running into numerical problems of overcom- 
pleteness. As a useful tool for the analysis of total energy data obtained from any 
first-principle calculation, we describe the scheme to solve the multi-dimensional vi- 
brational Schrodinger equation in the approximation of non-interacting anharmonic 
oscillators. The preliminary results for the lowest energy difference are presented 
for KNb0 3 . 
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